#include <stdio.h>
#include <math.h>
#include <string.h> 
#include <stdlib.h>
#include <time.h>
#define R 8.314
#define Comp_MAX 30  //define the numbers of the components
#define React_MAX 30 //define the numbers of reactions
#define Phase_MAX 5 //define the numbers of phase
void putfile();//double T,double a,double b,double c,double d,double e,double f);
void readfile();
void cal_abc();

void cal_delta();
void cal_i();
double cal_lnk298(int i); //i is the reaction number
double cal_dH_T(double T);
double cal_lnk(int i,double T);
double cal_lnk_ap(int i,double T);
void cal_eq_comp();
//void add_component();
//void add_reaction();
void cal_eq_comp_single();
void putfile_read();
void report();
void putfile_k(int tag,double k[10][3]);
void putfile_eq_comp(double T,double a,double b,double c,double d,double e,double f,int tag);
void calconstant();
void caldiagram();
void readfile_phasediag();
void calcp();


struct thermaldata
{
char name[20];
//double h,g,s,cp[5][2];//cp[i][0] for store the temperature, cp[i][1] for the heat capacity.
double cp[5][2];
int sta;//the state of the material.
double p[6];  // the parameter of heat capacity H298,G298,S298, a,b,c
}component[Comp_MAX+1];

/*char reactdef[React_MAX][2][50]={{4,1,1,4,4,3,3,6},
							{4,5,1,4,4,2,3,6},
							{4,2,3,6,4,5,1,4},
							{1,5,1,13,1,2,1,12},
							{1,2,1,12,1,5,1,13},
							{1,2,1,12,1,5,1,13},
							{1,2,1,12,1,5,1,13}};*/
struct reactiondata
{
char fullname[100];
int n_reactant;
int reactant_no[10];
int reactant_v[10];
int n_product;
int product_no[10];
int product_v[10];
double para[8];//8parameter, 0deltaH298,1deltaG298,2deltaS298,3deltaa,
//4deltab,5deltac,6deltaH0,7I .
}reaction[React_MAX],reaction_pd[Phase_MAX];

int comp_N,react_N,phase_N;
float n_CO_0=1,n_H2_0=1;//define the original CO and H2 mole number
float n_orig[React_MAX]={1,1,1,1,1};//define the original CO and H2 mole number
char errmsg[50];

void main()
{
	int code;
	readfile();
	cal_abc();
	cal_delta();
	do {    printf("\nWelcome to use the thermaldynamic calculating software\n\n	Please select the calculation type:\n	1 = The equibilium constant of reactions (reaction defined in reaction.txt)\n	2 = The equilibrium phase diagram calculation(reaction defined in PhaseDiagram.txt)\n	3=Caculate the Heat capcity (Cp) of component on the particular temperature.\n 0 = exit\nplease input the code:");
			scanf("%d",&code);
			switch(code)
			{
			case 0:
				return;
			case 1:
				calconstant();
				break;
			case 2:
				caldiagram();
				break;
			case 3:
				calcp();
				break;
			}
		}while(1);
}

void calconstant()
{
	int i,j;
double k[10][3];

//	printf("Please input the ratio of CO:H2, formart \"a:b\" eg. 1:1\n");
//	scanf("%f:%f",&n_CO_0,&n_H2_0);

	cal_i();
	putfile_read();
	for(i=0;i<react_N;i++)
	{
		for(j=0;j<10;j++)
		{
		k[j][0]=298.15+j*100;
		k[j][1]=cal_lnk(i,k[j][0]);
		k[j][2]=cal_lnk_ap(i,k[j][0]);
		};
//	putfile_k(i,k);
	}
	report();

 //	cal_eq_comp();
//	 cal_eq_comp_single();

//	printf("ok!");
//	getchar();
}
void caldiagram()
{
readfile_phasediag();
printf("\ncoming soon!\n");
}
void calcp()
{
int i,n;
char c;
float temp;
double cp;
//	printf("Please input the ratio of CO:H2, formart \"a:b\" eg. 1:1\n");
//	scanf("%f:%f",&n_CO_0,&n_H2_0);
do {    printf("The list of the component:\n");
		for(i=1;i<comp_N+1;i++)
		printf("%d\t%s\n",i,component[i].name);	
		do{
		printf("Please input the number of the component(0 for exit):");
		scanf("%d",&n);	
		if(n==0) return;		
		}while(n<1 || n>comp_N);
		
		printf("\nPlease input the temperature(K)(between 200k and 1000k):");
		scanf("%f",&temp);

		cp=component[n].p[3]+component[n].p[4]*temp+component[n].p[5]*temp*temp;
		printf("The heat capacity(Cp) of %s at %fk is %f\n continue(y/n):",component[n].name,temp,cp);
		scanf("\n%c",&c);
		if(c=='n') return;
		}while(1);



}
void cal_abc()
{

	int i,j,k;
	double p[3],x[3],y[3];

	for(i=1;i<comp_N+1;i++)
	{
		for(j=0;j<3;j++)
		{
		 x[j]=component[i].cp[2*j][0];
		 y[j]=component[i].cp[2*j][1];
		 if(y[j]==0)
		 {
			if(component[i].cp[2*j-1][1]!=0)
			 {
			 x[j]=component[i].cp[2*j-1][0];
			 y[j]=component[i].cp[2*j-1][1]; 
			 }
			else
			{
				strcpy(errmsg,"error");
				break;
			};
		 }
		}
		component[i].p[3]=0;
		component[i].p[3]=0;
		component[i].p[3]=0;
	 for(j=0;j<3;j++)
	 {
		 p[0]=1;
		 p[1]=0;
		 p[2]=1;
		 for(k=0;k<3;k++)
		 {
		 if(j!=k)
			{
			 p[0]=p[0]*x[k];
			 p[1]=p[1]+x[k];
			 p[2]=p[2]*(x[j]-x[k]);
			}
		 }
		component[i].p[3]=component[i].p[3]+p[0]*y[j]/p[2];
		component[i].p[4]=component[i].p[4]-p[1]*y[j]/p[2];
		component[i].p[5]=component[i].p[5]+y[j]/p[2];
	 }
	}
//	printf("cal_abc OK!\n");
}

void cal_delta()
{
	int i,j,k;
	for(i=0;i<react_N;i++)
	{
		for(j=0;j<6;j++)
		{   reaction[i].para[j]=0;
		for(k=0;k<reaction[i].n_reactant;k++)
			reaction[i].para[j]=reaction[i].para[j]
				-component[reaction[i].reactant_no[k]].p[j]*reaction[i].reactant_v[k];
		for(k=0;k<reaction[i].n_product;k++)
			reaction[i].para[j]=reaction[i].para[j]
				+component[reaction[i].product_no[k]].p[j]*reaction[i].product_v[k];
		}
	}
//	printf("cal_delta OK!\n");
}

void cal_i()
{
	int i;
	double dh;
	for(i=0;i<react_N;i++)
	{
		reaction[i].para[6]=reaction[i].para[0]*1000-reaction[i].para[3]*298.15
			-reaction[i].para[4]*298.15*298.15/2-reaction[i].para[5]*298.15*298.15*298.15/3;
		reaction[i].para[7]=cal_lnk298(i)+reaction[i].para[6]/R/298.15-reaction[i].para[3]/R*log(298.15)
			-reaction[i].para[4]/2/R*298.15-reaction[i].para[5]/6/R*298.15*298.15;
		dh=-reaction[i].para[6]/R/298.15+reaction[i].para[3]/R*log(298.15)
			+reaction[i].para[4]/2/R*298.15+reaction[i].para[5]/6/R*298.15*298.15+reaction[i].para[7];
	}
}

double cal_lnk298(int i)
{
	double result;

	result=-reaction[i].para[1]*1000/R/298.15;

return	result;//-reaction[i].para[1]*1000/R/298.15;
}

double cal_dH_T(int i,double T)
{
	double result;
	result=reaction[i].para[6]+reaction[i].para[3]*T+reaction[i].para[4]*T*T/2+reaction[i].para[5]*T*T*T/3;
	return result/1000.0;
}

double cal_lnk(int i,double T)
{
	double result;
	result=-reaction[i].para[6]/R/T+reaction[i].para[3]/R*log(T)+reaction[i].para[4]/2/R*T+reaction[i].para[5]/6/R*T*T+reaction[i].para[7];
	return result;
}

double cal_lnk_ap(int i,double T)
{
	double result,result1;
	result=-reaction[i].para[0]*1000/R*(1/T-1/298.15)+cal_lnk298(i);
	result1=exp(result);
	return result;
}

void cal_eq_comp()
{
	int i,tag=0;
	double T,epsilon_a,epsilon_b,K_a,K_b,yco2,yco,yh2o,yh2,N;
	N=n_CO_0+n_H2_0;
	for(i=0;i<11;i++){
	T=298.15+i*100;

	K_a=exp(cal_lnk(0,T));
	K_b=exp(cal_lnk(1,T));

	epsilon_a=n_CO_0*(1-1/(pow(K_a,0.25)+1))/4;
	epsilon_b=n_H2_0*(1-1/(pow(K_b,0.25)+1))/4;
	yco2=4*epsilon_a/N;
	yco=(n_CO_0-4*epsilon_a)/N;
	yh2o=4*epsilon_b/N;
	yh2=(n_H2_0-4*epsilon_b)/N;

	if(i==0)
		tag=1;
	else
		tag=0;
	putfile_eq_comp(T,epsilon_a,epsilon_b,yco2,yco,yh2o,yh2,tag);
	}
}

void cal_eq_comp_single()
{
	int i,j,tag=0;
	double T,epsilon,K,n_r,n_p;
	//N=n_CO_0+n_H2_0;
	for(i=0;i<react_N;i++)
	{
		for(j=0;j<11;j++){
		T=298.15+j*100;

		K=exp(cal_lnk(i,T));

		epsilon=n_orig[i]*(1-1/(pow(K,0.25)+1))/4;

		n_p=4*epsilon;
		n_r=n_orig[i]-4*epsilon;

		if(j==0)
			tag=1;
		else
			tag=0;
		putfile_eq_comp(T,epsilon,0,n_p,n_r,0,0,tag);
		}
	}
}
void readfile()
{
        FILE *fp;                
		int i,j,num;					//component 0 leave empty
		char *r_temp,*p_temp,*temp;
	 	fp=fopen("input.txt", "r");  
		   if(fp==NULL)            
             printf("File open error, input.txt doesn't exist ");

		   for(i=1;i<Comp_MAX+1;i++)
		   {
			  fscanf(fp,"%s\t%d\t%lf\t%lf\t%lf\t",
			   &component[i].name,&component[i].sta,&component[i].p[0],&component[i].p[1],&component[i].p[2]);
			  
			  for(j=0;j<4;j++)
				fscanf(fp,"%lf\t",&component[i].cp[j][1]);

				fscanf(fp,"%lf\n",&component[i].cp[4][1]);

				component[i].cp[0][0]=298.15;
				component[i].cp[1][0]=400.0;
				component[i].cp[2][0]=600.0;
				component[i].cp[3][0]=800.0;
				component[i].cp[4][0]=1000.0;
printf("debug component name %send\n",component[i].name);
				if(component[i].name==NULL)
				{printf("debugcomponent number %d\n",i);
					break;}
		   };
    comp_N=i-1;
printf("component number %d\n",comp_N);
i=fclose(fp);    

	 	fp=fopen("reaction.txt", "r");  
		   if(fp==NULL)            
             printf("File reaction.txt open error, input.txt doesn't exist ");

		   for(i=0;feof(fp)==0;i++)/*i<React_MAX;i++)*/
		   {
			  // fp=fp-1;
			   fscanf(fp,"%s\n",
				   &reaction[i].fullname);
	//	 printf ("%s\n",reaction[i].fullname);
		 temp=reaction[i].fullname;
/*		 		if(temp==NULL)
				   break;*/
			   r_temp = strtok (temp,"=");  
			 //  printf ("%s\n",r_temp);
			   p_temp = strtok(NULL,"=");
			 //  printf ("%s\n",p_temp);
			   j=0;
			   temp=strtok(r_temp,"+_");
			   while(temp!=NULL) {
				   num= atoi (temp);
				   if(j%2==0)
					   reaction[i].reactant_v[j/2]=num;
				   else
					   reaction[i].reactant_no[(j-1)/2]=num;
				   
				   temp= strtok(NULL,"+_");
				   j++;
			   };
			   reaction[i].n_reactant=j/2;
	 //  printf("reactant number %d",reaction[i].n_reactant);		   
			   j=0;
			   temp=strtok(p_temp,"+_");
			   while(temp!=NULL) {
				   num=atoi(temp);
				   if(j%2==0)
					   reaction[i].product_v[j/2]=num;
				   else
					   reaction[i].product_no[(j-1)/2]=num;
				   
				   temp= strtok(NULL,"+_");
				   j++;
			   };
			   reaction[i].n_product=j/2;
		//	      printf("product number %d\n",reaction[i].n_product);

				  if(i>React_MAX)
					  break;
		   };
   react_N=i ; 
   //printf("reaction number %d\n",react_N);
i=fclose(fp);         
//printf("readfile OK!\n");
}

void readfile_phasediag()
{
	 FILE *fp;                
		int i,j,num;					//component 0 leave empty
		char *r_temp,*p_temp,*temp;
	 	fp=fopen("PhaseDiagram.txt", "r");  
		   if(fp==NULL)            
             printf("File open error, input.txt doesn't exist ");

		    for(i=0;feof(fp)==0;i++)/*i<React_MAX;i++)*/
		   {
			  // fp=fp-1;
			   fscanf(fp,"%s\n",
				   &reaction_pd[i].fullname);
	//	 printf ("%s\n",reaction[i].fullname);
		 temp=reaction_pd[i].fullname;
/*		 		if(temp==NULL)
				   break;*/
			   r_temp = strtok (temp,"=");  
			 //  printf ("%s\n",r_temp);
			   p_temp = strtok(NULL,"=");
			 //  printf ("%s\n",p_temp);
			   j=0;
			   temp=strtok(r_temp,"+_");
			   while(temp!=NULL) {
				   num= atoi (temp);
				   if(j%2==0)
					   reaction_pd[i].reactant_v[j/2]=num;
				   else
					   reaction_pd[i].reactant_no[(j-1)/2]=num;
				   
				   temp= strtok(NULL,"+_");
				   j++;
			   };
			   reaction_pd[i].n_reactant=j/2;
	 //  printf("reactant number %d",reaction[i].n_reactant);		   
			   j=0;
			   temp=strtok(p_temp,"+_");
			   while(temp!=NULL) {
				   num=atoi(temp);
				   if(j%2==0)
					   reaction_pd[i].product_v[j/2]=num;
				   else
					   reaction_pd[i].product_no[(j-1)/2]=num;
				   
				   temp= strtok(NULL,"+_");
				   j++;
			   };
			   reaction_pd[i].n_product=j/2;
		//	      printf("product number %d\n",reaction[i].n_product);

				  if(i>Phase_MAX)
					  break;
		   };
   phase_N=i ; 
   //printf("reaction number %d\n",react_N);
i=fclose(fp);         
//printf("readfile OK!\n");
}

void putfile_read()
{ 
          FILE *fp;                
		  int i,j;
		  float co=n_CO_0,h2=n_H2_0;


	char filename[64];
	sprintf(filename, "Thermal(read).xls");

	fp=fopen(filename, "a+");  
		   if(fp==NULL)             
             puts("File open error"); 


     	for(i=1;i<comp_N+1;i++)
		{
		fprintf(fp,"%s\t%f\t%f\t%f\t",
		component[i].name,component[i].p[0],component[i].p[1],component[i].p[2]);
	
		for(j=0;j<5;j++)
			fprintf(fp,"%4.2f@%4.4f\t",component[i].cp[j][0],component[i].cp[j][1]);
		fprintf(fp,"%f\t%f\t%f\n",component[i].p[3],component[i].p[4],component[i].p[5]);

		}
		fprintf(fp,"deltaH298\tdeltaG298\tdeltaS298\tdelata\tdeltab\tdeltac\n");
		for(i=0;i<react_N;i++)
		fprintf(fp,"%f\t%f\t%f\t%f\t%f\t%f\n",reaction[i].para[0],
			reaction[i].para[1],reaction[i].para[2],reaction[i].para[3],
			reaction[i].para[4],reaction[i].para[5]);
i=fclose(fp);           
         if(i!=0)                  
            printf("File close error and %d",i); 
}

void putfile_k(int tag, double k[10][3])
{ 
          FILE *fp;                
		  int i;

	char filename[64];
	sprintf(filename, "Thermal(K).xls");

	fp=fopen(filename, "a+");  
		   if(fp==NULL)             
             puts("File open error"); 

	fprintf(fp,"reaction %d \t",tag);
	for(i=0;i<reaction[tag].n_reactant;i++)
	{
fprintf(fp,"%d %s",reaction[tag].reactant_v[i],component[reaction[tag].reactant_no[i]].name);
	if (i!=(reaction[tag].n_reactant-1))
		fprintf(fp,"+");
	}

	fprintf(fp,"->");

	for(i=0;i<reaction[tag].n_product;i++)
	{
fprintf(fp,"%d %s",reaction[tag].product_v[i],component[reaction[tag].product_no[i]].name);
	if (i!=reaction[tag].n_product-1)
fprintf(fp,"+");
	else
fprintf(fp,"\n");		
	}
     	for(i=0;i<10;i++)
		{
    fprintf(fp,"%f\t%f\t%f\n",k[i][0],k[i][1],k[i][2]);
		}

	i=fclose(fp);            
    if(i!=0)                  
     printf("File close error and %d",i); 
}

void putfile_eq_comp(double T,double a,double b,double c,double d,double e,double f,int tag)
{ 
          FILE *fp;                
		  int i;
		  float co=n_CO_0,h2=n_H2_0;


	char filename[64];
	sprintf(filename, "Thermal(CO_H2=%4.2f_%4.2f).xls", co,h2);

	fp=fopen(filename, "a+");  
		   if(fp==NULL)             
             puts("File open error");
		   if(tag==1)
				fprintf(fp,"T\tepsilonA\tepsilonB\tYCO2/n_p\tYCO/n_r\tYH2O\tYH2\n");
fprintf(fp,"%f\t%f\t%f\t%f\t%f\t%f\t%f\n",T,a,b,c,d,e,f);
     
i=fclose(fp);           
          if(i!=0)                  
           printf("File close error and %d",i); 
}      

void report()
{ 
          FILE *fp;                
		  int i,j;
		  float co=n_CO_0,h2=n_H2_0;
		  double T;
		  char filename[64],fullname[64];
		  struct tm *newtime; 
		  time_t t; 
			t = time(NULL); 
			newtime = localtime(&t); 
			strftime( filename, 64, "report%Y%m%d%H%M", newtime); 

	sprintf(fullname,"%s.txt",filename);
	
		fp=fopen(fullname, "r"); 
		   if(fp==NULL)             
		   {
			   fp=fopen(fullname, "w");
		   }
		   else
		   {
			   i=1;
			sprintf(fullname,"%s_%d.txt",filename,i);
		   while(fp=fopen(fullname, "r"))
		   {
		   i++;
		   sprintf(fullname,"%s_%d.txt",filename,i);
		   }
			fp=fopen(fullname, "w");
		   }
		   
				fprintf(fp,"This report is generated by Thermalequil.exe. This program can calculate the thermodynamic related data for reactions. More functions such as the equibilium conversion rate will present soon.");

				
		printf("\ndebug reaction%d\n\n",react_N);	
			for(i=0;i<react_N;i++)
			{	
			  fprintf(fp,"\n///////////////////////////////////////////////////////\n\n");
				fprintf(fp,"reaction %d \n",i);
				for(j=0;j<reaction[i].n_reactant;j++)
				{
				fprintf(fp,"%d %s",reaction[i].reactant_v[j],component[reaction[i].reactant_no[j]].name);
				if (j!=(reaction[i].n_reactant-1))
					fprintf(fp,"+");
				}

				fprintf(fp,"->");

				for(j=0;j<reaction[i].n_product;j++)
				{
				fprintf(fp,"%d %s",reaction[i].product_v[j],component[reaction[i].product_no[j]].name);
				if (j!=reaction[i].n_product-1)
				fprintf(fp,"+");
				else
				fprintf(fp,"\n");		
				}
				fprintf(fp,"delta H298 is % 4.5f(KJ/mol),\t",reaction[i].para[0]);
				fprintf(fp,"delta G298 is % 4.5f(KJ/mol),\t",reaction[i].para[1]);
				fprintf(fp,"delta S298 is % 4.5f(J/k/mol).\n",reaction[i].para[2]);
			
				if(reaction[i].para[0]>0)
					fprintf(fp,"This reaction is endothermic when the temperature is 298.5K\n");
				else
					fprintf(fp,"This reaction is exothermic when the temperature is 298.5K\n");
				fprintf(fp,"The data below show the effect of temperature on the reaction\n\n");
				fprintf(fp,"T(K)   \tdelta_H(KJ/mol)    \tlnK\n");
				for(j=0;j<20;j++)
				{
				T=198.15+j*50;
					fprintf(fp,"%4.2f\t%5.4f\t%5.4f\n",T,cal_dH_T(i,T),cal_lnk(i,T));
				};
			}
 
			printf("The report %s has been generated.\n You can view the result in the file.\nThank you for using.\n",fullname);
i=fclose(fp);           
          if(i!=0)                  
           printf("File close error and %d",i); 
}


